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Abstract 



We introduce calculus-based formulas for the continuous Euler and homotopy operators. The 
ID continuous homotopy operator automates integration by parts on the jet space. Its 3D 
generalization allows one to invert the total divergence operator. As a practical application, 
we show how the operators can be used to symbolically compute local conservation laws of 
nonlinear systems of partial differential equations in multi-dimensions. 

By analogy to the continuous case, we also present concrete formulas for the discrete Euler 
and homotopy operators. Essentially, the discrete homotopy operator carries out summation 
by parts. We use it to algorithmically invert the forward difference operator. We apply the 
discrete operator to compute fluxes of differential-difference equations in (1+1) dimensions. 

Our calculus-based approach allows for a straightforward implementation of the oper- 
ators in major computer algebra system, such as Mathematica and Maple. The symbolic 
algorithms for integration and summation by parts are illustrated with elementary examples. 
The algorithms to compute conservation laws are illustrated with nonlinear PDEs and their 
discretizations arising in fluid dynamics and mathematical physics. 

1 Introduction 

This chapter focuses on symbolic methods to compute polynomial conservation laws of par- 
tial differential equations (PDEs) in multi-dimensions and differential- difference equations 
(DDEs) (semi-discrete lattices). For the latter we treat only (1+1) dimensional systems 
where time is continuous and the spacial variable has been discretized. 

There are several strategies to compute conservation laws of PDEs. Some methods use a 
generating function [2], which requires the knowledge of key pieces of the Inverse Scattering 
Transform [1]. Other methods use Noether's theorem to get conservation laws from varia- 
tional symmetries. More algorithmic methods, some of which circumvent the existence of a 
variational principle [5, 6, 24, 33], require the solution of a determining system of ODEs or 
PDEs. Despite their power, only a few of these methods have been implement in computer 
algebra systems (CAS), such as Mathematica, Maple, and REDUCE. See [15, 33] for reviews. 

We advocate a more direct approach by building the candidate density as a linear com- 
bination (with constant coefficients) of terms that are uniform in rank (with respect to the 
scaling symmetry of the PDE). Although restricted to polynomial densities and fluxes, our 
method is entirely algorithmic and can be implemented in most CAS. We refer the reader to 
[15, 16] for details about an implementation in Mathematica, which can be downloaded from 
[18]. An implementation in Maple is also available [11]. 

Our earlier algorithm [15, 16] worked only for nonlinear PDEs in one spacial variable. In 
this chapter we present an algorithm that works for systems of PDEs in multi-dimensions 
that appear in fluid mechanics, elasticity, gas dynamics, general relativity, (magneto-)hydro- 
dynamics, etc. The new algorithm produces densities in which all divergences and equivalent 
terms have been removed. An additional advantage of our methods to compute densities and 
fluxes is that they can be applied to nonlinear DDEs [16, 19, 20]. 

During the development of our methods we came across tools from the calculus of varia- 
tions and differential geometry that deserve attention in their own right. These tools are the 
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variational derivative, the higher Euler operators, and the homotopy operator. 
To set the stage, we address a few issues arising in multivariate calculus: 

(i) To determine whether or not a vector field F is conservative, i.e. F = V/ for some scalar 
field /, one must verify that F is irrotational, that is V x F = 0. The field / can be computed 
via standard integrations [29, p. 518, 522]. 

(ii) To test if F is the curl of some vector field G, one must check that F is incompressible or 
divergence free, i.e. V • F = 0. The components of G result from solving a coupled system of 
first-order PDEs [29, p. 526]. 

(iii) To verify whether or not a scalar field / is the divergence of some vector function F, no 
theorem from vector calculus comes to the rescue. Furthermore, the computation of F such 
that / = V • F is a nontrivial matter. In single variable calculus, it boils down to computing 
the primitive F — J f dx. 

In multivariate calculus, all scalar fields /, including the components F t of vector fields 
F = (Fi, F 2 , F 3 ), are functions of the independent variables (x,y,z). In differential geome- 
try one addresses the above issues in much greater generality. There, the functions / and 
Fi can depend on arbitrary functions u(x,y, z),v(x,y, z), etc. and their mixed derivatives 
(up to a fixed order) with respect to the independent variables (x,y,z). Such functions are 
called differential functions [30]. As one might expect, carrying out the gradient-, curl-, or 
divergence-test requires advanced algebraic machinery. For example, to test whether or not 
/ = V • F requires the use of the variational derivative (Euler operator) in 3D. The actual 
computation of F requires integration by parts. That is where the homotopy operator and 
the variational complex come into play. 

In ID problems the continuous total homotopy operator 1 reduces the problem of symbolic 
integration by parts to an integration with respect to a single variable. In 2D and 3D, the 
homotopy operator allows one to invert the total divergence operator and, again, reduce 
the problem to a single integration. At the moment, no major CAS have reliable routines 
for integrating expressions involving unknown functions and their derivatives. As far as we 
know, no CAS offer a function to test if a differential function is a divergence. Routines to 
symbolically invert the total divergence are certainly lacking. 

The continuous homotopy operator is a universal, yet little known, tool that can be applied 
to many problems in which integration by parts (of arbitrary functions) in multi-variables 
plays a significant role. We refer the reader to [30, p. 374] for a history of the homotopy 
operator in the context of inverse problems of the calculus of variations. 

A major motivation for writing this chapter is to demystify the homotopy operators. 
Therefore, we purposely avoid differential forms and abstract concepts such as the variational 
bicomplex. Instead, we present down-to-earth calculus formulas for the homotopy operators 
which makes them readily implementable in major CAS. 

By analogy with the continuous case, we also present formulas for the discrete versions 
of the Euler and homotopy operators. The discrete homotopy operator is a powerful tool 
to invert the forward difference operator, whatever the application is. It circumvents the 
necessary summation (by parts) by applying a set of variational derivatives followed by a 
one-dimensional integration with respect to an auxiliary parameter. We use the homotopy 
operator to compute conserved fluxes of DDEs. Numerous examples of such DDEs are given 

1 Hence forth, homotopy operator instead of total homotopy operator. 
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in [31]. Beyond DDEs, the discrete homotopy operator has proven to be useful in the study 
of difference equations [22, 27]. To our knowledge, CAS offer no tools to invert the forward 
difference operator. Once fully implemented, our discrete homotopy operator will overcome 
the shortcomings. 

As shown in [22, 27], the parallelism between the continuous and discrete cases can be 
made rigorous and both theories can be formulated in terms of variational bicomplexes. 
To make our work accessible to as wide an audience as possible, we do not explicitly use 
the abstract framework. Aficionados of de Rham complexes may consult [7, 8, 9, 25] and 
[22, 27, 28]. The latter papers cover the discrete variational bicomplexes. 

2 Examples of Nonlinear PDEs 

We consider nonlinear systems of evolution equations in (3 + 1) dimensions, 
u t = G(u, 

Ux) Uy, U 2 , u.2xj **2y> **1zi "-xyi "-xzi ^-yzy ■ ■ ■ ), (1) 

where x = (x, y, z) are space variables and t is time. The vector u(x, y, z, t) has N components 
Ui. In the examples we denote the components of u by u, v, w, etc. Subscripts refer to partial 
derivatives. For brevity, we use u 2x instead of u xx , etc. and write G(u(™)) to indicate that 
the differential function G depends on derivatives up to order n of u with respect to x, y, and 
z. For simplicity, we assume that G does not explicitly depend on x and t. No restrictions 
are imposed on the number of components, the order, and the degree of nonlinearity of the 
variables in G. 

We will predominantly work with polynomial systems, although systems involving one 
transcendental nonlinearity can also be handled. If parameters are present in (1), they will 
be denoted by lower-case Greek letters. 

Example 1: The coupled Korteweg-de Vries (cKdV) equations [1], 

u t - Q(3uu x + 6vv x - f3u 3x = 0, v t + 3uv x + v 3x = 0, (2) 

where f3 is a nonzero parameter, describes interactions of two waves with different dispersion 
relations. System (2) is known in the literature as the Hirota-Satsuma system. It is completely 
integrable [1, 21] when f3 — \. 

Example 2 The sine-Gordon (SG) equation [10, 26], u-it — u 2x = sin-u, can be written as a 
system of evolution equations, 

u t = v, v t = u 2x + sin-u. (3) 

This system occurs in numerous areas of mathematics and physics, ranging from surfaces 
with constant mean curvature to superconductivity. 

Example 3: The (2+l)-dimensional shallow water wave (SWW) equations [12], 

Ui + (u-V)u + 2Q x u = - V(/i0) + \hV9, 

t + u-(V0) = 0, h t + V-(hu) = 0, (4) 
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describe waves in the ocean using layered models. Vectors u = u(x,y,t)i + v(x,y,t)j and 
Q = Q\t are the fluid and angular velocities, respectively, i, j, and k are unit vectors along 
the x, y, and z-axes. 9(x,y,t) is the horizontally varying potential temperature field, and 
h(x, y, t) is the layer depth. The dot (•) stands for Euclidean inner product and V = + J^j 
is the gradient operator. System (4) is written in components as 

u t + uu x + vu y — 2Qv + \h9 x + 9h x = 0, v t + uv x + vv y + 2Qu + \h6 y + 9h y = 0, 
9 t + u9 x + v9 y = 0, h t + hu x + uh x + hv y + vh y = 0. (5) 

3 Key Definitions-Continuous Case 

Definition: System (1) is said to be dilation invariant if it is invariant under a scaling 
(dilation) symmetry. 

Example: The cKdV system (2) is invariant under the scaling symmetry 

(x, t, u, v) -> (X^x, A" 3 t, \ 2 u, \ 2 v), (6) 
where A is an arbitrary scaling parameter. 

Definition: We define the weight, W , of a variable as the exponent of A that multiplies the 
variable. 

Example: We will always replace x by \^x. Thus, W(x) = —1 or W(d/dx) = 1. From (6), 
we have W(d/dt) = 3 and W(u) = W(v) = 2 for the cKdV equations. 

Definition: The rank of a monomial is defined as the total weight of the monomial. An 
expression (or equation) is uniform in rank if its monomial terms have equal rank. 

Example: Coincidentally, all monomials in both equations of (2) have rank 5. Thus, (2) is 
uniform in rank. The ranks of the equations in (1) may differ from each other. 

Conversely, requiring uniformity in rank for each equation in (1) allows one to compute the 
weights of the variables (and thus the scaling symmetry) with linear algebra. 
Example: For the cKdV equations (2), one has 

W{u) + W(d/dt) = 2W{u) + 1 = 2W(v) + 1 = W(u) + 3, 

W(v) + W(d/dt) = W{u) + W(v) + 1 = W(v) + 3, (7) 

which yields W(u) = W(v) = 2, W(d/dt) = 3, which leads to (6). 

Dilation symmetries, which are special Lie-point symmetries, are common to many nonlinear 
PDEs. However, non-uniform PDEs can be made uniform by extending the set of depen- 
dent variables with auxiliary parameters with appropriate weights. Upon completion of the 
computations one can set these parameters to one. 

Example: The sine-Gordon equation (3) is not uniform in rank unless we replace it by 

u t — v, v t = u 2x + asinu, a G 1R. (8) 
Using the Maclaurin series for the sin function, uniformity in rank requires 
W(u) + W(d/dt)=W(v), 

W(v) + W(d/dt) = W{u) +2 = W(a) +W{u) = W(a) + 3W(u) = W(a) +5W(u) = ■■■ . (9) 
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This forces us to set W(u) = 0. Then, W(a) = 2. By allowing the parameter a to scale, (8) 
becomes scaling invariant under the symmetry 

(x, t, u, v, a) — > (A _1 :r, A _1 t, X°u, X 1 v , A 2 a), (10) 

corresponding to W(d/dx) = W(d/dt) = l,W(u) = 0,W(v) = l,W(a) = 2. The first and 
second equations in (8) are uniform of ranks 1 and 2, respectively. 

Definition: System (1) is called multi-uniform in rank if it admits more than one dilation 
symmetry (which is not the result of adding auxiliary parameters with weights). 
Example: Uniformity in rank for the SWW equations (5) requires, after some algebra, that 

W(d/dt) = W(Q), W(d/dy) = W(d/dx) = 1, W{u) = W(v) = W(Q) - 1, 
W(6) = 2W(to)-W(h)-2, (11) 

where W(h) and IU(f2) remain free. The SWW system is thus multi-uniform. The symmetry 

(x, y, t, u, v, 9, h, n) -> (X^x, X^y, A~ 2 t, \u, Xv, X9, Xh, A 2 fi), (12) 

which is most useful for our computations later on, corresponds to W(d/dx) — W{d/dy) = 
l,W{d/dt) = 2,W(u)=W(v) = l,W(6) = l,W(h) = l, and W(Q) = 2. A second symmetry, 

(x, y, t, u, v, 9, h, Q) -> (X^x, X^y, X'H, Am, Xv, X 2 9, X°h, A 2 fi), (13) 

matches W (d / dx)=W (d / dy)=l , W(d/dt)=2, W(u)=W(v)=l, W(9)=2, W(h)=0, W(fi)=2. 

4 Conserved Densities and Fluxes of Nonlinear PDEs 

Definition: A scalar differential function p(u^) is a conserved density if there exists a vector 
differential function J(u^ m )), called the associated flux, such that 

D t p + DivJ = (14) 

is satisfied on the solutions of (1). Equation (14) is called a local 2 conservation law [30], 
and Div is called the total divergence 3 . Clearly, Div J = (D x , D y , D z ) ■ (Ji, J 2 , J3) = D x Ji + 
D y J 2 + D^J 3 . In the case of one spacial variable (x), (14) reduces to 

B t p + D x J = 0, (15) 

where both density p and flux J are scalar differential functions. In the ID case, 

where is the highest order term present in p. Upon replacement of ut,u tx , etc. from 
u t = G, one gets 

B tP = f t + p(uY[G], (17) 



2 We only compute densities and fluxes free of integral terms. 
3 Gradient, curl, and divergence are in rectangular coordinates. 
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where p(u)'[G] is the Frechet derivative of p in the direction of G. Similarly, 



a j m f) J 

D ^„( m , ) = _ + g_ M(i+1)i . (18) 

The generalization of (16) and (18) to multiple dependent variable is straightforward. For 
example, taking u = (u,v), 

Dtp( „<~»,„ M) = | + + (19) 

ft T 1711 r) T m2 ft T 

D x j(«<->y«*)) = m + ^m^<h (2Q) 



We will ignore densities and fluxes that explicitly depend on x and t. If G is polynomial then 

most, but not all, densities and fluxes are also polynomial. 

Example: The first four density-flux pairs for the cKdV equations (2) are 

p (1) = u, J (1) = -3/5w 2 + 3v 2 - [3u 2x , (any (3) (21) 

p( 2 ) = u 2 — 2v 2 , J (2) = -4/5w 3 + (3u 2 x - 2(3uu 2x + 2v 2 x - Avv 2x , (any (3) (22) 

p (3) = W, J (3) =3M 2 W + 2M 3 -M :E t; :E + M2xW + MW2z, (/3 = -1) (23) 



and 



p (4) = (l + (3)u 3 -3uv 2 -±(l + (3)u x 2 + 3v x 2 , (24) 
j( 4 ) = -1/^(1 + (3) u * + 9/3uV - \v A + 6/5(1 + /5)ww 2 - 3/5(1 + /?)u 2 u 2iE 
+3/3f 2 -u 2x — + + /^(l + P) u xU3x — §(3vu x v x + V2uv 2 x 
-6mw 2 s - 3u| x + 6v x v 3x ((3^-1). (25) 

The above densities are uniform in ranks 2,4 and 6. Both p^ and p^ 3 ^ are of rank 4. The 
corresponding fluxes are also uniform in rank with ranks 4, 6, and 8. In [15], we listed a few 
densities of rank > 8, which only exist when f3 — |. 

In general, if in (15) rankp = R then rank J = R + W(d/dt) — 1. All the terms in (14) 
are also uniform in rank. This comes as no surprise since the conservation law (14) holds on 
solutions of (1), hence it "inherits" the dilation symmetry of (1). 
Example: The first few densities [3, 13] for the sine-Gordon equation (8) are 

p {1) = 2acosu + v 2 + u x 2 , J (1) = -2u x v, (26) 

p {2) = 2u x v, J {2) = 2acosu-v 2 -u x 2 , (27) 

p^ = 6avu x cosw + v 3 u x + vu x — 8v x u 2x , (28) 

p( 4 ) = 2a 2 cos 2 u — 2a 2 sin 2 u + Aav 2 cos u + 20au x 2 cos u + v 4 + 6v 2 u x + u x 4 

-lQv x 2 - lQu 2 2x . (29) 

j( 3 ) and are not shown due to length. Again, all densities and fluxes are uniform in rank 
(before a is set equal to 1). 



6 



Example: The first few conserved densities and fluxes for the SWW equations (5) are 

p aw, J« = (j*), ^>=M, J'^(S), ^=M», J«=(^), (30) 
,«> = (u » + + k% ,«» = ( $ + ™^ + j p(5)= ^ _ v + 2W _ (31) 

(5) _ 1 / 12fiw6> - 4^0 + Quv x 9 + 2vv y 6 + u 2 9 y + v 2 9 y - h66 y + h y 6 2 N 



(32) 



6 v I2ttv9 + 4w x 6» - 6to^ - 2uu x 9 - u 2 9 x - v 2 9 x + h99 x - h x 9 

All densities and fluxes are multi-uniform in rank, which will substantially simplify the com- 
putation of the densities. Under either of the two scales, (12) or (13), rank(J) = rank(p) + I. 
With the exception of and j( 2 \ the ranks of the densities under scales (12) and (13) differ 
by one. The same holds for the fluxes. 



5 Tools from the Calculus of Variations 

In this section we introduce the variational derivative (Euler operator), the higher Euler op- 
erators from the calculus of variations, and the homotopy operator from homological algebra. 
These tools will be applied to the computation of densities and fluxes in Section 7. 



5.1 Continuous Variational Derivative (Euler Operator) 

Definition: A scalar differential function / is a divergence if and only if there exists a vector 
differential function F such that / = DivF. In ID, we say that a differential function / is 
exact if and only if there exists a scalar differential function F such that / = D X F. Obviously, 
F = D~ 1 (/) = J f dx is then the primitive (or integral) of /. 
Example: Consider 

f = 3u x v 2 sin-u — u x sinw — 6v v x cosw + 2u x U2 X cosw + 8^%, (33) 

which we encountered [3] while computing conservation laws for (8). / is exact. Indeed, upon 
integration by parts (by hand), one gets 

F = 4 v 2 + u 2 x cos u — 3v 2 cos u. (34) 

Most CAS, including Mathematica, Maple 5 and Reduce, fail this elementary integration! 
Example: Consider 

f = U X Vy - U 2x Vy - UyV X + U X yV X . (35) 

It is easy to verify that / = Div F with 

F = (uv y - u x v y , -uv x + u x v x ). (36) 

4 We do not use integrable to avoid confusion with complete integrability from soliton theory. 

5 Version 9.5 of Maple can integrate such expressions as a result of our interactions with the developers. 
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As far as we know, the leading CAS have no tools to compute F. Three questions arise: 

(i) Under what conditions for / does a closed form for F exist? 

(ii) If / is a divergence, what is it the divergence of? 

(hi) Avoiding integration by parts, how can one design a fast algorithm to compute F? 
To answer these questions we use the following tools from the calculus of variations: the 
variational derivative (Euler operator), its generalizations (higher Euler operators), and the 
homotopy operator. 

Definition: The variational derivative (Euler operator), 4°(L)' is defined [30, p. 246] by 

<>, - E(-d)^j. P?) 

where the sum is over all the unordered multi-indices J [30, p. 95]. For example, in the 
2D case the multi-indices corresponding to second-order derivatives can be identified with 
{2x, 2y, 2z, xy, xz, yz}. Obviously, (— D) 2x = D^, {—D) xy = D^D^, etc. For notational details 
see [30, p. 95, p. 108, p. 246]. 

With applications in mind, we give explicit formulas for the variational derivatives in ID, 2D, 
and 3D. For scalar component u they are 

4(1) = t(-V*) k ^- = #- - D * + Dl^- - Dl^- + • • • , (38) 
1 ' du kx du du x du 2x du 3x 



oo oo o 

,(o,o) _ Vr-D ) k *(-T) ) k » — 

'u(x,i/) ~ 1^ 2^ \ u x) { u y) ^ 



k y y 

= 7T - D *1T- -Vyi~ + d Itt- +D - D y7T- +D ?7T^ - D -7r ' ( 39 ) 

OU OU x OUy OU2 x OU xy OU2y OU 3x 



9 ~ 9 - 9 ^2 9 lT ^^ 9 lT ^2 9 ^3 9 



and 



oo oo oo 



4°£l) = E E E (-D.K-DyM-D,? 



k x =0 k y =0 k z =0 9u kxXkyykzZ 

d „ d „ d „ d „ d d „ 2 d 



D x — D„- D 2 — + Dt- + D- + D 



du x du x y du y z du z x du2 X y du2 V z du2 Z 
+ D X D W ^- + D X D Z ^- + D,D 2 ^- - Dl^- - • • • . (40) 

uu X y yju,y Z uu,3 x 

Note that u kxXkyy stands for u xx ... xyy ... y where x is repeated k x times and y is repeated k y 
times. Similar formulas hold for components v,w, etc. 

The first question is then answered by the following theorem [30, p. 248]. 

Theorem: A necessary and sufficient condition for a function / to be a divergence, i.e. 

there exists a F so that / = DivF, is that 4(x)(/) = 0- m other words, the Euler operator 

annihilates divergences, just as the divergence annihilates curls, and the curl annihilates 

gradients. 
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If, for example, u = (u,v) then both C^Ljf) and C^Ljf) must vanish identically. For the 
ID case, the theorem says that a differential function / is exact, i.e. there exists a F so that 
/ = D X F, if and only if C^ x) (f) = 0. 

Example: To test the exactness of / in (33) which involves just one independent variable 
x, we apply the zeroth Euler operator (37) to / for each component of u = (u, v) separately. 
For component u (of order 2), one computes 

4t,(/) = ^-D^ + D^ 



< X ^ J> du x du x x du 2x 

= 3u x v 2 cos u — v? x cos u + 6v v x sin u — 2u x u 2x sin u 

— Da;[3f 2 sin-u — 3u 2 sin-u + 2u 2x cos-u] + D 2 x \2u x cos-u] 
= 3u x v 2 cos u — ul cos u + 6v v x sin u — 2u x u 2x sin u 
— [3u x v 2 cos-u + 6vv x sin-u — 3-u^ cos-u — 6uu 2x smu 
—2u x u 2x sin u + 2u-$ x cos u] 
+ [— 2u 3x cos u — 6u x u 2x sin u + 2u 3x cos u] 
= 0. (41) 

Similarly, for component v (also of order 2) one readily verifies that >C^L(/) = 0. 
Example: As an example in 2D, one can readily verify that / = u x v y — u 2x v y — u y v x + 
u xy v x from (35) is a divergence. Applying (39) to / for each component of u = (u, v) gives 

5.2 Continuous Higher Euler Operators 

To compute F = Div _1 (/) or, in the ID case F = D~ x (/) = / f dx, we need higher-order 
versions of the variational derivative, called higher Euler operators. The general formulas are 
given in [30, p. 367]. With applications in mind, we restrict ourselves to the ID, 2D, and 3D 

cases. 

Definition: The higher Euler operators in ID (with variable x) are 

(4 2 ) 

where (^fj is the binomial coefficient. Note that the higher Euler operator for % = matches 
the variational derivative in (38). The explicit formulas for the first three higher Euler 
operators (for component u and variable x) are 



^» = ^- 3D ^ + 6D ^- 10D ^ + ''-- (44) 

£ «"> = <L - 4D -i + ioD2 ^ - 2oD ^ + ■ - < 45 > 
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Definition: The higher Euler operators in 2D (with variables x, y) are given by 



oo oo 



^Si = E Ell (-D^-(-D^^--^. (46) 
k x =i x ky=i y \ Lx / \ L y / vu.k xX k y y 



Note that the higher Euler operator for i x — i y — matches the variational derivative in (39). 
The first higher Euler operators (for component u and variables x and y) are 

r^o) - A 9D — p _ I qp 2 ® | 2D D - I D 2 — (AT) 

£ ^ = du~ y ~ 2By du7y- 3D ^ + 2DKD ^ + d *^t ■ ■ ■ ' (48) 



4(^) = J— 2 ^7T- 2D,-^-+3D 2 -^-+4D :c D s/ -^-+-- - , (49) 

du xy du 2xy du x2y du 3xy du 2x2y 



£ S) = ^ 3D *^ 2D ^^+ 6D ^+ 3D ^ • ( 5 °) 

OU 2xy OU 3xy OU 2x2 y OU^y OU 2x3 y 

Definition: The higher Euler operators in 3D (with variables x, y, z) are 

c$i$=t £ £ f t ¥* !, y tN )(-D.)'--(-D ! ,)'=«-(-D I )'-- " 



"^.^AW WW 1 " 1 1 a«^MM' (5l) 

The higher Euler operator for i x = i y = i z = matches the variational derivative given in (40). 



5.3 Continuous Homotopy Operator 

We now discuss the homotopy operator which will allow us to reduce the computation of 
F = Div _1 (/) (or in the ID case, F = D~ 1 (f) = / f dx) to a single integral with respect to an 
auxiliary variable denoted by A (not to be confused with A in Section 3). Hence, the homotopy 
operator circumvents integration by parts and reduces the inversion of the total divergence 
operator, Div, to a problem of single-variable calculus. The homotopy operator is given in 
explicit form, which makes it easier to implement in CAS. To keep matters transparent, we 
present the formulas of the homotopy operator in ID, 2D, and 3D. 

Definition: The homotopy operator in ID (with variable x) [30, p. 372] is 

««(.)(/)= rXX(/)[Au]^, (52) 
Jo j=l A 

where the integrand I Uj (/) is given by 

oo 

M/)=E D i (53) 

i=0 

The integrand involves the ID higher Euler operators in (42). In (52), N is the number of 
dependent variables and /„.(/) [Au] means that in I Uj (f) one replaces u(x) — > \u(x), u x (x) — > 
Xu x (x), etc. 
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Given an exact function /, the question how to compute F = D x 1 (f) = J f dx is then 
answered by the following theorem [30, p. 372]. 
Theorem: For an exact function /, one has F = Tt u (x)(f)- 

Thus, in the ID case, applying the homotopy operator (52) allows one to bypass integration 
by parts. A clever argument why the homotopy operator actually works is given in [6, p. 
582]. As an experiment, one can start from a function F, compute / = D X F, subsequently 
compute F = Ti u {x)(f), an d finally verify that F — F is a constant. 

Example: Using (33), we show how the homotopy operator (52) is applied. For a system 
with N = 2 components, (ui,u 2 ) = (u,v), the homotopy operator formulas are 

W»(x)(/) = jf 1 (/«(/) [Aii] + W)[Au]) y , (54) 

with 

OO CO 

W) = E D x(«4w(/)) and W) = E D «(«4w 1) (/))- (ss) 

i=0 i=0 

These sums have only finitely many non-zero terms. For example, the sum in I u (f) terminates 
at p — 1 where p is the order of u. Take, for example, f = 3u x v 2 sinu—u x sin u — 6 v v x cos u + 
2u x u 2x cosm + 8v x v 2x . First, we compute 

W) = «£? ( L)(/) + D x (u£« } (/)) 

u- 2uD x I - | + B x I u 



OU x \OU 2x J \ OU 2x J 

3uv 2 sin u — uu x sin u + 2u x cos u. (56) 



Next, 



W) = <l)(/) + D*K(i)(/)) 



cto x \ov 2x ) \ av 2x ) 

= -Qv 2 cosu + 8vl. (57) 

Formula (54) reduces to an integral with respect to A : 

f = n u(x) (f) = J\w)[\u] + i v (f)[\u}) y 

= / (3X 2 uv 2 sin(Xu) — X 2 uu 2 sin(Xu) + 2Xu 2 ,cos(Xu) — 6Xv 2 cos(Xu) + 8Xv 2 ^jdX 

J 

= Av 2 + u 2 x cos u — 3v 2 cos u. (58) 
We now turn to the problem of inverting the Div operator using the homotopy operator. 
Definition: We define the homotopy operator in 2D (variables x, y) through its two compo- 
nents (Ti.u{x ^u(L y)(f))- The x-component of the operator is given by 

= L £4?(/)[Au]y, (59) 

JO =1 A 
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with Ijf(f) given by 

_l+i x 

i x =0 i y =0 \ 

Analogously, the y-component is given by 



4f(f) = EEJ^-j ^ («i4^(/)) • ( 6 °) 



«2[L,y)(/) = / n E^(/)[Au]y, (61) 

JU ,=1 A 



with 

1 + / 

=0i u =0 ' 



These integrands involve the 2D higher Euler operators in (46). 

After verification that / is a divergence, the question how to compute F = (Fi,F 2 ) = 
Div _1 (/) is then answered by the following theorem. 

Theorem: If F is a divergence, then F = (F 1 ,F 2 ) = Div"^/) = (n^ y) (f),H^ y) (f)). 

The superscript (x) in 1iS x \f) reminds us that we are computing the ^-component of F. As a 
test, one can start from any vector F and compute / = DivF. Next, compute F = (F 1; F 2 ) = 
(K%, y) (f),H%, y) (f)) and, finally, verify that F - F is divergence free. 

Example: Using (35), we show how the application of the 2D homotopy operator leads to 
(36), up to a divergence free vector. Consider / = u x v y — u 2x v y — u y v x + u xy v x , which is easily 
verified to be a divergence. In order to determine Div _1 (/), we calculate 



df OT 3/ 3/ \ , / 0/ \ 1/ 3/ 



Similarly, 



Hence, 



U \du x 2 ^~ >x du 2x ^ >y du xy ) + ^ x \ du 2x ) + 2^ y \du xy ) 

UVy + ~UyV X - U X Vy + ^UV X y . (63) 



= ^4S)(/) = «J£ = "V + V- (64) 
= wil s) (/) = /(4 x) (/)[Au]+/W(/)[Au]) y 

= J A ^UUj, + -Wj/Ur - + ~UV xy - U y V + U xy vj d\ 



111111 

-MVy + ~UyV X - ~U x Vy + ~UV X y - -UyV + ~U X yV. (65) 
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Without showing the details, one computes in an analogous fashion 

F2 = WS[L >l ,)(/)=jf 1 (/S' ) (/)[Au]+/^(/)[Au]) y 

= J A (-uv x - ^uv 2x + ^u x v x ^j + A (u x v - u 2x v) dX 

11111 

= -^uv x - -uv 2x + -u x v x + -u x v - -u 2x v. (66) 

One can readily verify that the resulting vector 

y — ( Fl \ = ( ^ uv y + ^ u y Vx ~ * UxV v + ^ UVx y ~ 5 V + ^ u *y v \ (67) 

\F 2 ) \ -\uv x - \uv 2x + \u x v x + \u x v - \u 2x v J 

differs from F = (uv y — u x v y , —uv x + u x v x ) by a divergence free vector. 
The generalization of the homotopy operator to 3D is straightforward. 
Definition: The homotopy operator in 3D (with variables x, y, z) is 

B y anal °gy with ( 59 )' the component is 

<U)(/) = / E^ } (/)[Au]y, (68) 

with, 



00 00 00 



iif(f) = E E E ( tttt^m ) D * D * D " (^<^r } (/)) • (69) 



i x =o i v =o i z =o \ x + % * + *y + 1 



The 7/ and ^-operators are defined analogously. The integrands involve the 3D higher Euler 
operators in (51). By analogy with the 2D case the following theorem holds. 

Theorem: Given a divergence / one has F = Div" 1 (/) = {H%^ z) {j\ H%^ z) {j\ H%^ z) {j)). 



6 Removing Divergences and Equivalent Terms 

We present an algorithm to remove divergences and equivalent terms in order to make our 
computation of the densities simpler. 

Definition: Two scalar differential functions, and f^ 2 \ are equivalent if and only if they 
differ by the divergence of some vector V, i.e. /« ~ if and only if - = Div V. 
Obviously, if a scalar expression is equivalent to zero, then it is a divergence. 
Example: Functions = uu 2x and = —u 2 x are equivalent because /W — = uu 2x + 
u 2 x = D x (uu x ). Using (38), note that v\ = C^(uu 2x ) = 2u 2x and v 2 = £u( x )(~ u l) = 2u 2x are 
equal. Also, / = u± x = D x {u^ x ) is a divergence and, as expected, t>3 = C^) x Ju4x) = 0. 

Example: In the 2D case, = (u x — u 2x )v y and = (u y — uxy)v x are equivalent since 
- f {2) = u x v y - u 2x v y - u y v x + u xy v x = Div (uv y - u x Vy, -uv x + u X V x ) . Using (39), note 

that Vl = C^^fW) = V 2 = C^l^fW) = (~V X y - V XX y, ~U X y + U XX y) . 

To remove divergences and equivalent terms we use the following algorithm. 
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Algorithm: Remove-Divergences- And-Equivalent-Terms (7?.) 



/* Given is list 1Z of monomial differential functions */ 
/* Initialize two new lists S, B */ 

/* Find first member of S */ 
for each term U E 1Z 

do Vi <- £$ X) (U) 
ifv^O 

then 5 <- {ti} 

break 
else discard t; and V; 
/* Find remaining members of 5 */ 
for each term tj G 1Z \ {t 1: t 2) • ■ • , U} 

ifv^O 

then if Vj ^ Span(£>) 

then S<-SU{tj} 

else discard i,- and Vj 

return S 

/* List S is free of divergences and equivalent terms */ 

Example: Using the above algorithm, we remove divergences and equivalent terms in 
TZ = {u 3 ,u 2 v,uv 2 ,v 3 ,u 2 x ,u x v x ,vl,uu 2x ,u 2x v,uv 2x ,vv 2x ,u 4x ,v 4x }. Since vi = C^(u 3 ) = 
(3m 2 , 0) ^ (0,0) we have S = {h} = {u 3 } and B = {vi} = {(3m 2 , 0)}. The first for loop 
is halted and the second for loop starts. Next, v 2 = C^^(u 2 v) = (2uv,u 2 ) ^ (0, 0). We ver- 
ify that Vi and v 2 are independent and update the sets resulting in S = {t l5 1 2 } = {u 3 ,u 2 v} 
and B = {vi, v 2 } = {(3m 2 , 0), (2uv, u 2 )}. 

Proceeding in a similar fashion, since the first seven terms are indeed independent, we 
have S = {h,t 2 , ■■■ ,t 7 } and B = {vi, v 2 , • • • , v 7 } = {(3m 2 , 0), (2uv,u 2 ), • • • , (0, -2v 2x )}. 
For t$ = uu 2x we compute v 8 = ^^(mm^) = {2u 2x , 0) and verify that v 8 = — v 5 . So, v 8 G 
Span(i3) and t s and v 8 are discarded (i.e. not added to the respective sets). For similar reasons, 
t g , t w , and tn as well as v 9 , Vi , and Vn are discarded. The terms t 12 = u ix and t ls = v± x are 
discarded because v 12 = v 13 = (0,0). So, 1Z is replaced by S = {u 3 , u 2 v , uv 2 , v 3 , m 2 , u x v x , v 2 } 
which is free of divergences and equivalent terms. 



7 Application: Conservation Laws of Nonlinear PDEs 

As an application of the Euler and homotopy operators we show how to compute conserved 
densities and fluxes for the three PDEs in Section 2. The first PDE illustrates the ID case 
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(one independent variable), but it involves two dependent variables u(x) and v(x). The second 
PDE (again in ID) has a transcendental nonlinearity which complicates the computation of 
conserved densities and fluxes [13]. A third example illustrates the algorithm for a 2D case. 

To compute conservation laws D t p + Div J = of systems of nonlinear PDEs, we use 
a direct approach. First, we build the candidate density p as a linear combination (with 
constant coefficients q) of terms which are uniform in rank (with respect to the scaling 
symmetry of the PDE). It is of paramount importance that the candidate density is free of 
divergences and equivalent terms. If such terms were present, their coefficients could not 
be determined because such terms can be moved into the flux, J. To construct the shortest 
density, we will use the algorithm of Section 6. 

Second, we evaluate D t p on solutions of the PDE, thus removing all time derivatives 
from the problem. The resulting expression (called E) must be a divergence (of the as yet 
unknown flux). Thus, we set ^)^a(E) = 0. Setting the coefficients of like terms to zero leads 
to a linear system for the undetermined coefficients q. In the most difficult case, such systems 
are parameterized by the constant parameters appearing in the given PDE. If so, a careful 
analysis of the eliminant (and solution branching) must be carried out. For each branch, the 
solution of the linear system is substituted into p and E. 

Third, since E = Div J we use the homotopy operator 7i u ( x ) to compute J = Div _1 (i£). 
The computations are carried out with our Mathematica packages [18] . 



7.1 Conservation Laws for the Coupled KdV Equations 

In (21) through (24) we gave the first four density-flux pairs. As an example, we will compute 
density p^ and associated flux J^. 

Recall that the weights for the cKdV equations are W(d/dx) = 1 and W(u) = W(v) = 2. 
The parameter (3 has no weight. Hence, p^ has rank 6. The algorithm has three steps: 

Step 1: Construct the form of the density 

Start from V = {u, v}, i.e. the list of dependent variables with weight. Construct M. = 
{u 3 ,v 3 ,u 2 v,uv 2 ,u 2 ,v 2 ,uv,u,v,l}, which contains all monomials of selected rank 6 or less 
(without derivatives). Next, for each monomial in Ai, introduce the correct number of x- 
derivatives so that each term has rank 6. For example, 

d 2 u 2 d 2 v 2 d 2 {uv\ 

-q^2~ = 2u l + 2uu 2x , -^j = 1v\ + 2vv 2x , ^2 = U2xV + 2UxVx + UV2x > 

d 4 u d 4 v d 6 l n ,_ n . 

W =Ui - W =Vi - ^ = °- (70) 

Ignore the highest-order terms (typically the last terms) in each of the right hand sides of 
(70). Augment M. with the remaining terms, after stripping off numerical factors, to get 
1Z = {u 3 ,u 2 v,uv 2 ,v 3 ,u 2 ,u x v x ,v 2 .,U2 X v}, where the 8 terms are listed by increasing order. 

Note that keeping all terms in (70) would have resulted in the list 71 (with 13 terms) given 
in the example at the end of Section 6. As shown, the algorithm would reduce 1Z to 7 terms. 

Use the algorithm of Section 6, to replace 1Z by S = {u 3 , u 2 v, uv 2 , v 3 , u 2 , u x v x , v 2 }. Linearly 
combine the terms in S with constant coefficients to get the shortest candidate density: 

p = civ 3 + c 2 u 2 v + c 3 uv 2 + c 4 v 3 + c 5 ul + c 6 u x v x + c 7 v 2 x . (71) 
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Step 2: Determine the constants q 

Compute 

dp . dp dp dp dp 

E — D t p = — + p (u) F = T~Ut + ~^—u tx + —v t + -^—v tx 
dt du du x dv dv x 

= (3ci« 2 + 2c 2 uv + c 3 v 2 )u t + (2c 5 u x + c e v x )u tx + (c 2 u 2 + 2c 3 uv + 3c 4 v 2 )v t 

+(c 6 u x + 2c 7 v x )v tx . (72) 

Replace u t ,v t ,u tx and v tx from (2) to obtain 

E = (3c 4 u 2 + 2c 2 uv + c 3 v 2 )(6/3uu x - 6vv x + (3u 3x ) + (2c 5 u x + c % v x )(Q(3uu x - 6vv x + (3u 3x ) x 
-(c 2 u 2 + 2c 3 uv + 3c 4 v 2 )(3uv x + v 3x ) - (c 6 u x + 2c 7 v x )(3uv x + v 3x ) x . (73) 

Since E = D t p = — D X J, the expression E must be exact. Therefore, apply the variational 
derivative (38) and require that C^ X ^(E) = and cf^(E) = 0. Group like terms and set 
their coefficients equal to zero to obtain the following (parameterized) linear system for the 
unknown coefficients c\ through c 7 : 

(3 + 4/3)c 2 = 0, 3 Cl + (l+/3)c 3 = 0, 4c 2 + 3c 4 = 0, (1 + (3)c 3 - 6c 5 = 0, 

/3(ci + 2c 5 ) = 0, /3c 2 -c 6 = 0, (l + /3)c 6 = 0, c 4 + c 6 = 0, (74) 

2(1 + /3)c 2 - 3(1 + 2/3)c 6 = 0, 2c 2 - (1 + 6/3)c 6 = 0, /3c 3 - 6c 5 - c 7 = 0, c 3 + c 7 = 0. 

Investigate the eliminant of the system. In this example, there exists a solution for any 
P 7^ — 1. Set ci = 1 and obtain 

ci = 1, c 2 = c 4 = c 6 = 0, c 3 = -j^p, c 5 = -|, c 7 = y|^. (75) 

Substitute the solution into (71) and multiply by 1 + (3 to get 

p = (1 + /3)w 3 - 3m; 2 - |(1 + (3)u x 2 + 3v x 2 , (76) 

which is p (4) in (24). 

Step 3: Compute the flux J 

Compute the flux corresponding to p in (76). Substitute (75) into (73), reverse the sign and 
multiply by 1 + /3, to get 

E = 18/3(1 + (3)v?u x - 18(3u 2 vv x - 18(3uu x v 2 + 18v 3 v x - 6/3(1 + /3)w 3 - 6/3(1 + (3)uu x u 2x 
+3/3(1 + f3)u 2 u 3x - 3(3v 2 u 3x - 6v x v 4x - /3(1 + (3)u x u 4x + Quvu 3x + 6(/3 - 2)u x v 2 x 
+6(1 + (3)u x vv 2x - \8uv x v 2x (77) 

Apply (52) and (53) to (77) to obtain 

J = -1(3(1 + (3)u A + 9(3u 2 v 2 - ft; 4 + 6/3(1 + (3)uu 2 x - 3/3(1 + (3)u 2 u 2x 
+3(3v 2 u 2x - \(3(l + (3)u\ x + (3(1 + (3)u x u 3x - 6/3vu x v x 
+12uv 2 - 6uvv 2x - 3v 2 x + 6v x v 3x , (78) 

which is in (25). 

The cKdV equations (2) are completely integrable if (3 — \ and admit conserved densities 
at every even rank. 
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7.2 Conservation Laws for the sine-Gordon Equation 

Recall that the weights for the sine- Gordon equation (3) are W{£) = 1, W{u) = 0, W(v) = 1, 
and = 2. The first few (of infinitely many) densities and fluxes were given in (26) 

through (29). We show how to compute densities p^ and p^ 2 \ both of rank 2, and their 
associated fluxes and J^. 

In contrast to the previous example, the candidate density will no longer have constant 
undetermined coefficients q but functional coefficients hi{u) which depend on the transcen- 
dental variable u with weight zero [3]. To avoid having to solve PDEs, we tacitly assume that 
there is only one dependent variable with weight zero. 

Step 1: Construct the form of the density 

Augment the list of dependent variables with a (with non-zero weight) and replace u by u x 
(since W{u) = 0). Hence, V = {a,u x ,v}. Compute 1Z = {a,v 2 ,v 2 ,u 2x ,u x v,u x } and remove 
divergences and equivalent terms to get S = {a,v 2 ,ul,u x v}. The candidate density 

p = ahiiu) + h 2 {u)v 2 + h 3 (u)u x 2 + h±{u)u x v , (79) 

with undetermined functional coefficients hi{u). 
Step 2: Determine the functions hi(u) 
Compute 

E = D tP = p'(u) [F] = ^u t + ^u tx + % t 

ot ou ou x ov 

= (ah[+v 2 h' 2 +u 2 ,h' 3 +u x vh' 4 )v + (2u x h 3 +vh i )v x +(2vh2+u x h4)(a sin(-u)+-u 2 a;)- (80) 

where h! i means Since E = D t p = —D X J, the expression E must be exact. Therefore, 

require that C^) \\{E) = and C^) x s{E) = 0. Set the coefficients of like terms equal to zero to 
get a mixed linear system of algebraic and ODEs: 

h 2 (u) -h 3 (u) = 0, ti 2 (u)=0, h' 3 (u) = 0, h' 4 (u)=0, h' 2 \u)=0, (81) 
h'l(u) = 0, 2h' 2 (u) - h' 3 (u) = 0, 2h' 2 \u) - h 3 (u) = 0, (82) 
h[(u) + 2h 2 (u) sinw = 0, h'[(u) + 2h' 2 (u) sinu + 2h 2 (u) cosu = 0. (83) 

Solve the system [3] and substitute the solution 

hi(u) = 2ci cos-u + c 3 , h 2 (u) — h 3 (u) — c±, h^u) — c 2 , (84) 

(with arbitrary constants q) into (79) to obtain 

p — Ci (2a cos u + v 2 + u x 2 ) + c 2 u x v + c 3 a. (85) 

Step 3: Compute the flux J 

Compute the flux corresponding to p in (85). Substitute (84) into (80), to get 

E = C\(2u 2x v + 2u x v x ) + c 2 (vv x + u x u 2x + au x sin-u). (86) 
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Since E — D t p = — D x J, one must integrate f — —E. Applying (55) yields I u (f) = —2c\U x v — 
c 2 (u x + au sin u) and I v (f) = —2ciu x v — c 2 t> 2 . Use formula (54) to obtain 



dX 

T 



J = Wu(x)(/)= f 1 (W)[\U\ + I v (f)[\u\) 
J 

= — J (Aci\u x v + C2(\u x + au sin(Aw) + Xv 2 )) dX 

= — Ci(2u x v) — c 2 (^v 2 + -u x — a cos u^j . (87) 

Finally, split density (85) and flux (87) into independent pieces (for c\ and c 2 ): 

p (1) = 2acosu + v 2 + u x 2 and J (1) = -2u x v, (88) 
p^ = u x v and = —-v 2 — ]-u 2 x + acosu. (89) 

For E in (86), J in (87) can easily be computed by hand [3]. However, the computation of 
fluxes corresponding to densities of ranks > 2 is cumbersome and requires integration with 
the homotopy operator. 



7.3 Conservation Laws for the Shallow Water Wave Equations 

In contrast to the previous two examples, as far as we know, (5) is not completely integrable. 
One cannot expect infinitely many conserved densities and fluxes (of different ranks). 

The first few densities and fluxes were given in (30). We show how to compute densities 
p^\ p( 3 \ p( 4 \ and p^\ which are of rank 3 under the following (choice for the) weights 

W(d/dx) = W(d/dy) = 1, W(u) = W(v) = 1, W{9) = 1, W{h) = 1, W(Q) = 2. (90) 

We will also compute the associated fluxes J( 3 \ J^, and J^. 

The fact that (5) is multi-uniform is advantageous. Indeed, one can use the invariance 
of (5) under one scale to construct the terms of p, and, subsequently, use additional scale(s) 
to split p into smaller densities. This "divide and conquer" strategy drastically reduces the 
complexity of the computations. 

Step 1: Construct the form of the density 

Start from V = {u, v, 9, h, Q}, i.e. the list of variables and parameters with weights. Use (90) 
to construct M. = {flu, flv, ■ ■ ■ , u 3 , v 3 , • • • , u 2 v, uv 2 , ■ ■ ■ ,u 2 , v 2 , uv, ■ ■ ■ ,u,v,9, h}, which has 
38 monomials of rank 3 or less (without derivatives). 

The terms of rank 3 in M. are left alone. To adjust the rank, differentiate each monomial 
of rank 2 in Ai with respect to x ignoring the highest-order term. For example, in ^j- = 2uu x , 
the term can be ignored since it is a total derivative. The terms u x v and —uv x are equivalent 
since = u x v + uv x . Keep u x v. Likewise, differentiate each monomial of rank 2 in M. with 
respect to y and ignore the highest-order term. 

Produce the remaining terms for rank 3 by differentiating the monomials of rank 1 in 
Ai with respect to x twice, or y twice, or once with respect to x and y. Again ignore the 
highest-order terms. Augment the set Ai with the derivative terms of rank 3 to get 1Z = 
{flu, flv, ■ ■ ■ , uv 2 , u x v, u x 9, u x h, ■ ■ ■ , u y v, u y 9, ■ ■ ■ , 9 y h} which has 36 terms. 
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Instead of applying the algorithm of Section 6 to TZ, use the "divide and conquer" strategy 
to split 7Z into sublists of terms of equal rank under the (general) weights 

W(d/dt) = W(n), W(d/dy) = W(d/dx) = 1, W(u) = W(v) = W(tt) - 1, 
W(9) = 2W(n)-W(h)-2, (91) 

where W(Q) and W(h) are arbitrary. Use (91), to compute the rank of each monomial in 1Z 
and gather terms of like rank in separate lists. 

Apply the algorithm from Section 6 to each TZi to get the list Si. Coinci dentally, in this 
example TZi = Si for all %. Linearly combine the monomials in each list Si with coefficients to 
get the shortest candidate densities p«. In Table 1, we list the 10 candidate densities and the 
final densities with their ranks. 



i 


Rank 


Candidate pi 


Final pi 


Final Jj 


1 


6W(Q)-3W(h)-6 


Cl e 3 








2 


3W(h) 


cih 3 








3 


hW{Vt)-2W{h)-h 


c 1 u9 2 +c 2 v9 2 








4 


W{Q)+2W{h)-l 


C\uh 2 +C2vh 2 








5 


AW{9)-W{h)-A 


c 1 u 2 9+c 2 uv9+c 3 v 2 9+c 4 9 2 h 


9 2 h 


/ uh9 2 \ 
{ vh9 2 ) 


6 


2W(Q)+W(h)-2 


CiU 2 h+C2Uvh-\-c 3 v 2 h-\-c i 9h 2 


u 2 h+v 2 h+9h 2 


f u 3 h+uv 2 h+2uh 2 9 \ 
\ v 3 h+u 2 vh+2vh 2 9 ) 


7 


3W(Q)-W(h)-2 


ciQ9+c 2 Uy9-\-C3Vy9+C4U x 9+C5V x 9 


2tt9-u y 9+v x 9 


J 7 in (102) 


8 


w(n)+w(h) 


CiQh^2Uyh-U: 3 Vyh-\-C4U x h-\-c 5 v x h 


Qh 


/ fluh \ 
{ Qvh ) 


9 


2W(tt)-l 


C\ Qu+C2^lv+C 3 U y V+C^yh+C^UxV+C^ x h 








10 


3W{Q)-3 


CiU 3 +C2U 2 V+C 3 UV 2 +C4V 3 +C5u9h+CQv9h 









Table 1: Candidate densities for the SWW equations. 
Step 2: Determine the constants q 

For each of the densities pi in Table 1 compute E; L = Y) t Pi an d use (5) to remove all time 
derivatives. For example, proceeding with p 7 , 



*,=„-(») [F] = + fi t% + pv„ + 3%„ + me, 

ou x ou y ov x OVy 09 

= —C/${uu x + vUy — 2Vtv + \h9 x + 9h x ) x — c 2 9{uu x + vu y — 2Qv + \h9 x + 9h x ) y 
—c§9{uv x + vv y + 2Qu + \h9 y + 9h y ) x — c 3 9(uv x + vv y + 2Qu + \h9 y + 9h y ) y 
-(ciQ + c 2 u y + c 3 v y + c 4 u x + c 5 v x )(u9 x + v9 y ). (92) 



dpi 



Require that C^] y) {E 7 ) = £^"> y) (E 7 ) = £^> y) (E 7 ) = £ ( ^> y) (E 7 ) = 0, where, for example, 
^u(x\) * s gi ven m (39). Gather like terms. Equate their coefficients to zero to obtain 

ci + 2c 2 = 0, c 3 = c 4 = 0, ci - 2c 5 = 0, c 2 + c 5 = 0. (93) 



•(0,0) 



•(0,0) 



.(0,0) 
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Set ci = 2. Substitute the solution 

ci = 2, c 2 = -1, c 3 = c 4 = 0, c 5 = 1. (94) 

into p 7 to obtain p 7 = 2Vt9 — u y 9 + 1^, which corresponds to p^ in (30). 

Proceed in a similar way with the remaining nine candidate densities to obtain the results 
given in the third column of Table 1 . 

Step 3: Compute the flux J 

Compute the flux corresponding to all pi ^ in Table 1. For example, continuing with p 7 , 
substitute (94) into (92) to get 

E 7 = -6{u x v x + uv 2x + v x v y + vv xy + 2Qu x + \9 x h y - u x u y - uu xy - u y v y - u 2y v 

+2Qv y — \6 y h x ) — (2Qu9 x + 2Qv9 y — uu y 6 x — u y v9 y + uv x 6 x + vv x 6 y ). (95) 

Apply the 2D homotopy operator in (59)-(61) to E 7 = — Div J 7 . So, compute 
iP(Eh) = uC^ y) (E 7 ) + D, (uC^ y) (E 7 )) + ±D W (uC^ y) {E 7 )) 

u( dEr 2D ( dE A d f^ 7 ^ D ^ dE 7 \ 1 / d£ 7 \ 
V^Ux H^J V \ du xyJJ \ du 2x ) 2 y \ du xy ) 

= -uv x 9 - 2Qu9 - -u 2 6 y + uu y Q. (96) 

2 



Similarly, compute 



i[ x \e 7 ) = - VVy e-^ v 2 e y -uv x e, (97) 

/^(Ey) = -^9 2 hy-2Qu9 + uu y e -uv x 9, (98) 
(£ 7 ) = Ue y h. (99) 



Next, compute 



j^(u) = -n { :l v) {E 7 ) 



- [ (l ( :\E7)[M + 4 x \E 7 )[\u] + lf\E 7 )[\u] + I^(E 7 )[\u}) ^ 

£ (i\Qu9 + X 2 (zuv x 9 + ^u 2 9 y - 2uu y d + vv y 9 + ^v 2 9 y + ^9 2 h y - ^y^)) 

2 11111 

2fiw6» - —uu y 9 + uua.0 + -w„0 + -u 2 9 y + -v 2 ^ - -/i00„ + -h y 9 2 . (100) 

3 3 6 6 6 6 



In analogous fashion, compute 



2fiw0 + -w x - vu v 9 - -uu x 9 - -u 2 9 x - -v 2 9 x + -h99 x - -h x 9 2 . (101) 
3 3 6 6 6 6 
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Hence, 



if nviue - Auu y e + &uv x e + 2w y e + u 2 e y + v 2 e y - hee y + h y e 2 \ 

7 ~ 6 V l2 ^v9 + Avv x 6 - 6vu y 9 - 2uu x - u 2 6 x - v 2 9 x + hOO x - h x 9 2 ) ' 

(102) 

which matches J^ 5 ^ in (32). 

Proceed in a similar way with the remaining nonzero densities to obtain the fluxes given 
in the last column of Table 1 . 

System (5) has conserved densities of the form 

p = hf(6) and p=(v x -u y + 2n)g(6), (103) 

where / and g are arbitrary functions. Our algorithm can only find polynomial / and g. A 
comprehensive study of all conservation laws of (5) is beyond the scope of this chapter. 



8 Examples of Nonlinear DDEs 

We consider nonlinear systems of DDEs of the form 

u n = G(..., u n _i, u„, u„+i, ...), (104) 

where u n and G are vector-valued functions with iV components. The integer n corresponds 
to discretization in space 6 ; the dot denotes differentiation with respect to the continuous time 
(t). For simplicity, we write G(u n ), although G depends on u„ and a finite number of its 
forward and backward shifts. We assume that G is polynomial with constant coefficients. 
No restrictions are imposed on the forward or backward shifts or the degree of nonlinearity 
in G. In the examples we denote the components of u n by u n ,v n , etc. If present, parameters 
are denoted by lower-case Greek letters. We use the following two DDEs to illustrate the 
theorems and algorithms. 

Example 4: The Kac-van Moerbeke (KvM) lattice [23], 

il n = U n (u n+ i - Un-i), (105) 

arises in the study of Langmuir oscillations in plasmas, population dynamics, etc. 
Example 5: The Toda lattice [32] in polynomial form [16], 

Un = v n -i - v n , v n = v n (u n -u n+1 ), (106) 

models vibrations of masses in a lattice with an exponential interaction force. 

9 Dilation Invariance and Uniformity in Rank for DDEs 

The definitions for the discrete case are analogous to the continuous case. For brevity, we 
use Example 5 to illustrate the definitions and concepts. 

6 We only consider DDEs with one discrete variable. 
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As shown in Table 3, the Toda lattice (106) is invariant under the scaling symmetry 



(t, u n , v n ) -> (A H, Xu n , X 2 v n ). 



(107) 



Definition: The weight W of a variable equals the exponent of the scaling parameter A. 
[16, 17]. Weights of dependent variables are nonnegative and rational. We tacitly assume 
that weights are independent of n. For example, W(u n -i) = W(u n ) = W(u n +i), etc. 

Example: Since t is replaced by | we have W(^) = W(D t ) = 1. From (107) we have 
W{u n ) = 1 and W(v n ) = 2. 

Definition: The rank of a monomial equals the total weight of the monomial. An expression 
(or equation) is uniform in rank if all its monomial terms have equal rank. 

Example: The three terms in the first equation in (106) have rank 2; all terms in the second 
equation have rank 3. Each equation is uniform in rank. 

Conversely, requiring uniformity in rank in (106) yields W(u n ) + 1 — W(v n ), and W{v n ) + 
1 — W(u n ) +W(v n ). Hence, W(u n ) — 1, W(v n ) — 2. So, the scaling symmetry can be computed 
with linear algebra. 

Many integrable nonlinear DDEs are scaling invariant. If not, they can be made so by 
extending the set of dependent variables with parameters with weights. 

10 Conserved Densities and Fluxes of Nonlinear DDEs 

By analogy with D x and D" 1 , we define the following operators acting on monomials m n in 

Uni etc. 

Definition: D is the upshift operator (also known as the forward- or right-shift operator) 
D m n = m n+1 . Its inverse, D -1 , is the down-shift operator (or backward- or left-shift operator), 
D _1 m n = m n _i. The identity operator is denoted by I, lm n = m n and A = D — I, is the 
forward difference operator. So, Am„ = (D — I) m n = m n+1 — m n . 

Definition: A conservation law of (104), 



which holds on solutions of (104), links a conserved density p n to a flux J n . Densities and 
fluxes depend on u n as well as forward and backward shifts of u n . 

Definition: Compositions of D and D _1 define an equivalence relation (=) on monomials. 
All shifted monomials are equivalent. 

Example: For example, -u n _it> n+ i = u n v n+2 = u n+ iv n+3 = u n+ 2V n+i . Factors in a monomial 
in u n and its shifts are ordered by u n+ j -< u n+ k if j < k. 

Definition: The main representative of an equivalence class is the monomial with u n in the 
first position [16, 17]. 

Example: The main representative in class {• • • , u n -2U n , Un-iUn+i, u n u n+ 2, ■ ■ ■} is u n u n+ 2 
(not u n _ 2 u n ). 

For monomials involving u n ,v n ,w n , etc. and their shifts, we lexicographically order the 
variables, i.e. u n -< v n -< w n , etc. For example, u n v n+2 (not ■u n _2f n ) is the main representative 

Of {• • • , U n - 2 V n , U n - 1 V n+1 , U n V n+2 , U n+1 V n+3 , ••■}. 



B tPn + AJ n = 0, 



(108) 
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To stress the analogy between PDEs and DDEs, we put the defining equations next to each 
other in Table 2. 





Continuous Case (PDE) 


Semi-discrete Case (DDE) 


Evolution Equation 


u t = G(u, u x ,u y , ...,u 2x , •••) 


u n =G(..., u n _i, u„, u n+ i, ...) 


Conservation Law 


B tP + V • J = 


p n + A J n = 



Table 2: Defining equations for conservation laws for PDEs and DDEs. 



Table 3 shows the KvM and Toda lattices with their scaling invariance (and weights) and a 
few conserved densities. Notice that the conservation law "inherits" the scaling symmetry of 
the DDE. Indeed, observe that all p n in Table 3 are uniform in rank, be it of different ranks. 





Kac-van Moerbeke Lattice 


Toda Lattice 


Lattice 


u n = u n (u n+ i - U n -i) 


U n = V n -! - V n , V n = V n (u n - U n+ i) 


Scaling 


(t,u n ) -> (A _1 t, \u n ) 


(t,u n ,v n ) -> (\~H, \u n ,\ 2 v n ) 


Weights 


W(D t ) = 1, W(u n ) = 1 


W(D t ) = 1, W(u n ) = 1, W(v n ) = 2 


Densities 


Pn l) = U n, P ( n } = U nU n +l 

Pf 1 = \ U n + U nU n +l{u n + U n+l +U n+2 ) 


P [ n ] = U n, P^ = \ul + V n 
P { n = \ U l + U n (v n -1 + Vn) 



Table 3: Examples of nonlinear DDEs with weights and densities. 



11 Discrete Euler and Homotopy Operators 

11.1 Discrete Variational Derivative (Euler Operator) 

Given is a scalar function f n in discrete variables u n ,v n , . . . and their forward and backward 
shifts. The goal is to find the scalar function F n so that f n = AF n = F n+1 — F n . We illustrate 
the computations with the following example: 

fn = -U n U n+ iV n -vl + U n+ iU n+2 V n+ i + V 2 n+l + U n+3 V n+2 ~ U n+1 V n . (109) 

By hand, one readily computes 

F n = v 2 n + u n u n+l v n + u n+1 v n + u n+2 v n+1 . (110) 
Below we will address the questions: 

(i) Under what conditions for /„ does F n exist in closed form? 

(ii) How can one compute F n = A _1 (/„) ? 

(iii) Can one compute F n = A _1 (/„) in an analogous way as in the continuous case? 
Expression /„ is called exact if it is a total difference, i.e. there exists a F n so that 

/„ = AF n . With respect to the existence of F n in closed form, the following exactness 
criterion is well-known and frequently used [4, 20]. 

Theorem: A necessary and sufficient condition for a function f n , with positive shifts, to be 
exact is that £g](/„) = 0. 
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is the variational derivative (discrete Euler operator of order zero) [4] defined by 

42 = ^-(E D ") = ^-(I + D " + D- 2 + D-3 + ■ ■ ■ ). (Ill) 

A proof of the theorem is given in e.g. [20]. In practice, the series in (111) terminates at 
the highest shift occurring in the expression the operator is applied to. To verify that an 
expression E(u n - q , • • • , u n , • • • , u n+p ) involving negative shifts is a total difference, one must 
first remove the negative shifts by replacing E n by E n = D q E n . 

Example: We return to (109), /„ = -u n u n+l v n -vl + u n+1 u n+2 v n+1 + w 2 +1 + u n+3 v n+2 - 
u n+iV n . We first test if f n is exact (i.e., the total difference of some F n to be computed 
later). We apply the discrete zeroth Euler operator to f n for each component of u n = (u n , v n ) 
separately. For component u n (with maximum shift 3) one readily verifies that 

4°2(/n) = ^- (I + D- 1 + D- 2 + D- 3 ) (fn) = 0. (112) 

Similarly, for component v n (with maximum shift 2) one checks that C^(f n ) = 0. 



11.2 Discrete Higher Euler and Homotopy Operators 

To compute F n , we need higher-order versions of the discrete variational derivative. They 
are called discrete higher Euler operators £W i n analogy with the continuous case [30]. 
In Table 4, we have put the continuous and discrete higher Euler operators side by side. Note 
that the discrete higher Euler operator for % = is the discrete variational derivative. The 
first three higher Euler operators for component u n from Table 4 are 

£« = A( D -i + 2D- 2 + 3D- 3 + 4D- 4 + ---), (113) 

£g) = (p~ 2 + 3D~ 3 + 6D~ 4 + 10D" 5 + •••). (114) 

ou n 

£g) = (p- ? > + 4D" 4 + 10D~ 5 + 20D" 6 + •••). (115) 



Similar formulae hold for C$. 





Continuous Case 


Discrete Case 


Zeroth Euler Operator 


rW v^oo ( T\ \k d 

^ u ( x )—2^k=0\ ^x) du (k) 


r(0)—\^co r\-k d 
L>u n -Lk=0 U d u n+k 


Higher Euler Operators 


/*(») _v^oo (k\( p> \k-i d 
^u(x)—2^k=i\ij\ u x) g u (k) 




Homotopy Operators 




W«.(/)=/oEf = i/« JlB (/)[Au B ]^ 


Integrand Operator 







Table 4: Continuous and discrete Euler and homotopy operators in ID side by side. 
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Also in Table 4, we put the formulae for the discrete homotopy operator T~t Un and the con- 
tinuous homotopy operator side by side. The integrand I Ujn (f) of the homotopy operator 
involves the discrete higher Euler operators. As is the continuous case, N is the number of 
dependent variables and I Uj n (f)[\u n ] means that after I Ujn (f) is applied one replaces u n by 
Au n , u n+ i by Au n+ i, etc. To compute F n , one can use the following theorem [19, 22, 27]. 
Theorem: Given an exact function f n , one can compute F n = A~ 1 (f n ) from F n = Tt Un (f n ). 

Thus, the homotopy operator reduces the inversion of A (and summation by parts) to a set of 
differentiations and shifts followed by a single integral with respect to an auxiliary parameter 
A. We present a simplified version [19] of the homotopy operator given in [22, 27], where 
the problem is dealt with in greater generality and where the proofs are given in context of 
variational complexes. 

For a system with components, {u\, n ,U2, n ) = (u n ,v n ), the discrete homotopy operator 
from Table 4 is 

n Un (f) = £ (AJ/MAun] + /J/)N) ^, (lie) 

with 

oo oo 

IuM) = E A * Mt +1) (/)) and /„„(/) =^A ! (y n C^\f)) . 

i=o i=o x (117) 

Example: We return to (109). Using (117), 

I Un (fn) = U n C^ n (fn) + A (u n £% (fn)) + A 2 (u n C% (fn)) 

= u n ^- (D- 1 + 2D" 2 + 3D- 3 ) (/„) + A [u n ^- (D- 2 + 3D- 3 ) (fn)} 

+ A 2 (Unl-B-Vn)) 

= 2u n u n+ iv n + u n+1 v n + u n+2 v n+1 , (118) 

and 

/»„(/,,) = «n£ii'(/») + A (»„££>(/»)) 

= ^( D " ,+2D " 2 )(« +A ("»l: D " 2 («) 

= u n u n+1 v n + 2v 2 n + u n+1 v n + u n+2 v n+1 . (119) 
The homotopy operator (116) thus leads to an integral with respect to A : 

Fn = £ (I Un (fn)[Xu n ] + I Vn (fn)[Xu n }) y 

= J (2Av 2 + 3\ 2 u n u n+1 v n + 2\u n+1 v n + 2\u n+2 v n +i) d\ (120) 
= v 2 n + u n u n+1 v n + u n+1 v n + u n+2 v n+1 , (121) 
which agrees with (110), previously computed by hand. 
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12 Application: Conservation Laws of Nonlinear DDEs 



In [14, 20], different algorithms are presented to compute ffuxes of nonlinear DDEs. In this 
section we show how to compute ffuxes with the discrete homotopy operator. For clarity, 
we compute a conservation law for Example 5 in Section 8. The computations are carried 
out with our Mathematica packages [18]. The completely integrable Toda lattice (106) has 
infinitely many conserved densities and fluxes. As an example, we compute density pffl 
(of rank 3) and corresponding flux J^f 1 (of rank 4). In this example, G = (G± : G 2 ) = 
( v n -i—v n , v n (u n —u n+ i) ). Assuming that the weights W(u n ) = 1 and W(v n ) — 2 are computed 
and the rank of the density is selected (say,i? = 3), our algorithm works as follows: 

Step 1: Construct the form of the density 

Start from V = {u n ,v n }, i.e. the list of dependent variables with weight. List all monomials 
in u n and v n of rank 3 or less: M. = {u n 3 , u n 2 , u n v n , u n , v n }. 

Next, for each monomial in Ai, introduce the correct number of t-derivatives so that each 
term has rank 3. Using (106), compute 

d°M n 3 _ 3 d°u n v n _ 
dt° ~ Un ' dt° ~ UnVn ' 

7] = zu n u n = ZU n V n —i Zu n v n , — — = V n = U n V n U n +iV n , 

dt dt 

d 2 U n diL n d(f n _i - V n ) 



dt 2 dt dt 



ti„-lt>n-l - u n v n - X - u n v n + U n+1 V n . (122) 



Augment M. with the terms from the right hand sides of (122) to get 
TZ = {M n 3 ,M n w n _i,M„t; n ,M„_it; n _i,M n+ it; n }. 

Identify members belonging to the same equivalence classes and replace them by their 
main representatives. For example, u n v n -\ = u n+ iv n , so the latter is replaced by u n v n -\. 
Hence, replace TZ by S = {u n 3 , u n v n -i, u n v n }, which has the building blocks of the density. 
Linearly combine the monomials in S with coefficients Cj to get the candidate density: 

p n = Ci Un + C 2 U n V n _i + C 3 U n V n . (123) 

Step 2: Determine the coefficients 

Require that (108) holds. Compute D t p n . Use (106) to remove u n and v n and their shifts. 
Thus, 

E n = D t p n = (3ci - c 2 )ulv n -i + (c 3 - 3ci)m^„ + (c 3 - c 2 )v n ^ x v n 

+ c 2 u n _ x u n v n _ x + c 2 vl_ 1 - c 3 u n u n+1 v n - c 3 vl. (124) 

To remove the negative shift n — 1, compute E n = DE n Apply to E n , yielding 

4°2(K) = ^-(I + D^ + D- 2 )^) 

= 2(3ci - c 2 )u n v n -i + 2(c 3 - 3ci)u n v n + (c 2 - c 3 )w n _it; n _i 

+ (c 2 - c 3 )u n+1 v n . (125) 
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Next, apply £g) to E n , yielding 

= ^-(I + D- 1 )^) 

= (3ci - c 2 )u 2 n+1 + (c 3 - c 2 )v n+ i + (c 2 - c 3 )u n u n+1 + 2(c 2 - c 3 )v„ 

+ (c 3 - 3ci)m^ + (c 3 - C2)u„_i. (126) 

Both (125) and (126) must vanish identically Solve the linear system 

3ci - c 2 = 0, c 3 - 3ci = 0, c 2 - c 3 = 0. (127) 

Set ci = | and substitute the solution c\ = |, c 2 = c 3 = 1, into (123) 

Pn = ^«n 3 +Wn(v„-1 + V n ). (128) 

Step 3: Compute the flux 

In view of (108), one must compute J n = — A~ 1 (E n ). Substitute c± — |,c 2 = c 3 = 1 into 
(124). Then, = D£ n = w n w n+1 t; n + v%- u n+1 u n+2 v n+1 - v% +1 . 
Apply (117) to — E n to obtain 

I Un {-E n ) = 2u n u n+1 v n , I Vn (-E n ) = u n u n+ iv n + 2v 2 n . (129) 

Application of homotopy operator (116) yields 

Jn = f 1 (I Un (-E n )[\u n ] + I Vn (-E n )[Xu n }) ^ 
Jo A 

= / (3\ 2 u n u n+1 v n + 2Xv^) dX 
Jo 

= u n u n+l v n + v 2 n . (130) 
After a backward shift, J n = D _1 (J n ), we obtain J n . With (128), the final result is then 

Pn = ^ul + U n (v n - 1 +V n ), J n = U n -iU n V n -i + V 2 n _ x . (131) 

The above density corresponds to in Table 3. 



13 Conclusion 

Based on the concept of scaling invariance and using tools of the calculus of variations, 
we presented algorithms to symbolically compute conservation laws of nonlinear polynomial 
and transcendental systems of PDEs in multi-spacial dimensions and DDEs in one discrete 
variable. We covered the symbolic computation of densities and fluxes. 

The continuous homotopy operator turned out to be a powerful, algorithmic tool to com- 
pute fluxes explicitly. Indeed, the homotopy operator handles integration by parts in multi- 
variables which allowed us to invert the total divergence operator. Likewise, the discrete 
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homotopy operator handles summation by parts and inverts the forward difference operator. 
In both cases, the problem is reduced to an explicit integral from ID calculus. 

Homotopy operators have a wide range of applications in the study of PDEs, DDEs, fully 
discretized lattices, and beyond. We extracted the continuous and discrete Euler and ho- 
motopy operators from pure mathematics, introduced them into applied mathematics, and 
therefore make them readily applicable to computational problems. We purposely avoided 
differential forms [30] and abstract concepts from differential geometry and homological al- 
gebra. 

Our down-to-earth approach might appeal to scientists who prefer not to juggle exte- 
rior products and Lie derivatives. Our calculus-based formulas for the Euler and homotopy 
operators can be readily implemented in major CAS. 
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